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The  sensitivity  of  a  rarehed-to-transitional  flow  to  the  fidelity  of  the  chemical  reaction  model  is 
investigated  for  a  new  molecular  dynamics/quasiclassical  trajectory  (MD/QCT) -derived  model  and 
compared  with  the  widely  used  total  collision  energy  (TCE)  model  of  Bird.  For  hypervelocity 
collisions  that  occur  in  the  space  environment,  it  is  not  clear,  a  priori,  that  the  TCE  model  will 
provide  reasonable  results  for  the  required  high  energy  range  and,  particularly,  if  strong  favoring  of 
the  reaction  among  different  forms  of  reactant  energy  occurs.  In  fact,  in  previous  work,  the  TCE 
model,  using  available  Arrhenius  parameters,  has  been  found,  for  these  flow  conditions,  to  give 
unphysical  probabilities.  A  chemical  reaction  model,  suitable  for  use  in  the  direct  simulation  Monte 
Carlo  (DSMC)  method,  is  developed  to  simulate  the  hypervelocity  collisions  of  0(3F) 
+  HC1(  *S+)  -f  OH(  2I1)  +  C1(2P),  an  example  of  an  important  reaction  in  high-altitude 
atmospheric-jet  interactions.  The  model  utilizes  the  MD/QCT  method  with  a  new  benchmark  triplet 
A"  surface.  Since  the  modeling  of  chemical  reactions  in  DSMC  simulations  requires  the  use  of  a 
reaction  probability,  the  adequacy  of  the  overall  collision  cross  section,  usually  modeled  by  the 
variable  hard  sphere  (VHS)  model,  is  also  considered.  To  obtain  an  accurate  collision  cross  section, 
the  approach  of  Tokumasu  and  Matsumoto  was  used  in  the  MD/QCT  method  with  the 
aforementioned  potential  energy  surface.  Energy  transfer  between  the  target  HCI  translational  and 
internal  energy  modes  was  investigated  and  it  was  found  that  the  variation  of  the  inelastic  cross 
section  has  a  negligible  effect  on  the  transport  cross  section.  Therefore,  a  MD/QCT  VHS  equivalent 
collision  cross  section  was  obtained  and  along  with  the  MD/QCT  reaction  cross  sections  were 
utilized  in  the  full  DSMC  calculation  of  the  flow  field.  It  was  found  that  for  a  low  enthalpy  reaction, 
in  hypervelocity  collisions,  the  TCE  model  with  accurate  Arrhenius  rates  appears  to  agree  well  with 
the  rigorous  MD/QCT  calculations  which  shows  that  the  reaction  does  not  exhibit  strong 
favoring.  ©  2007  American  Institute  of  Physics.  [DOI:  10.1063/1.2717692] 


I.  INTRODUCTION 

Treatment  of  chemical  reactions  in  rarefied,  nonequilib¬ 
rium  flows  has  been  the  subject  of  significant  research.  In 
some  cases,  proper  treatment  of  reactions  is  needed  to  accu¬ 
rately  predict  overall  shock  thickness,  standoff,  and  heat 
transfer.1  In  other  cases,  such  as  when  reactants  are  minor 
species,  the  overall  flow  field  is  not  affected  by  details  of  the 
reaction  model,  but  a  particular  radiation  signature  may  be 
strongly  affected.  A  case  of  the  latter  situation  is  the  subject 
of  the  present  study,  and  allows  for  a  detailed  examination  of 
several  issues  important  to  chemical  reaction  models  for  rar¬ 
efied,  high-speed  flows.1'4 

The  most  widely  used  method  for  calculating  reaction 
probabilities  in  direct  simulation  Monte  Carlo  (DSMC)  (Ref. 
5)  is  the  total  collision  energy  (TCE)  model.5  6  Collision 
pairs  are  sampled  according  to  the  usual  method,  and  each 
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pair  is  tested  for  reaction.  The  reaction  probability  is  an  as¬ 
sumed  function  of  the  total  collision  energy  of  the  reactant 
molecules  that,  under  equilibrium  conditions,  will  reproduce 
experimental  reaction  rates  KfT)  in  the  modified  Arrhenius 
form, 

Kf=ATn  e\p(-  Ea/kT),  (1) 

where  A  is  the  pre-exponential  factor,  Ea  is  the  activation 
energy,  and  n  is  an  additional  temperature  dependence.  In 
addition,  the  TCE  model  requires  an  input  parameter  defin¬ 
ing  the  number  of  degrees  of  freedom  of  the  internal  energy 
modes  of  the  reactants  that  contribute  to  the  reaction  prob¬ 
ability. 

There  are  two  major  concerns  with  the  TCE  chemistry 
model.  First,  while  the  assumed  form  of  the  probability  re¬ 
produces  measured  rates  under  equilibrium  conditions,  it 
may  not  accurately  describe  the  reaction  under  highly  non¬ 
equilibrium  conditions,  especially  if  the  reaction  is  activated 
primarily  via  energy  in  a  particular  mode  (i.e.,  favoring) 
rather  than  via  the  total  energy.  Second,  while  literature  rate 
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expressions  (such  as  those  to  be  discussed  in  this  work)  are 
fit  to  data  and  are  valid  only  for  a  prescribed  range  of  tem¬ 
peratures,  high  velocity  flows  involve  collision  energies  that 
correspond  to  a  much  higher  temperature  range.  Extrapolat¬ 
ing  Arrhenius  parameters  which  are  only  valid  up  to  3000  K 
into  the  10  000  K  range  and  beyond  should  be  avoided  if  at 
all  possible.  For  instance,  extrapolating  the  expression  of 
Mahmud7  (defined  below)  beyond  7000  K  leads  to  reaction 
rates  that  were  found  in  earlier  work4  to  be  higher  than  the 
total  collision  rate — a  clearly  unphysical  situation.  In  gen¬ 
eral,  the  TCE  model  provides  reasonable  results  if  it  is  used 
with  sound  Arrhenius  parameter  values  for  the  relevant  en¬ 
ergy  range  and  if  there  is  no  strong  favoring  of  the  reaction 
among  different  forms  of  reactant  energy.  While  one  cannot 
criticize  the  use  of  this  approach  if  no  better  information  is 
available,  it  certainly  motivates  a  search  for  alternatives. 

The  proposed  approach  to  be  discussed  in  this  work  ad¬ 
dresses  both  concerns  mentioned  above.  The  molecular 
dynamics/quasiclassical  trajectory  (MD/QCT)  method  of 
computing  reaction  cross  sections  can  be  expected  to  provide 
reliable  results  over  the  energy  range  of  interest,  so  long  as 
the  potential  energy  surface  (PES)  used  is  adequate  and 
quantum  mechanical  effects  such  as  tunneling  are  not  impor¬ 
tant.  The  emphasis  of  this  paper  is  to  evaluate  the  adequacy 
of  the  relatively  simple  TCE  model  for  the  O  +  HCl— >  OH 
+  C1  reaction,  important  in  high-altitude  atmospheric -jet  in¬ 
teractions.  Chemical  reaction  probabilities  obtained  by  the 
TCE  method  may  then  be  tested  with  this  more  fundamental 
approach  involving  new  calculations  using  the  MD/QCT 
method  and  the  recent,  benchmark  triplet  A"  PES  for  O 

o 

+  HC1.  Since  the  modeling  of  chemical  reactions  in  DSMC 
simulations  requires  the  use  of  a  reaction  probability,  the 
adequacy  of  the  overall  collision  cross  section,  usually  mod¬ 
eled  by  the  variable  hard  sphere  (VHS)  model,  must  also  be 
considered.  To  obtain  an  accurate  collision  cross  section,  the 
approach  of  Tokumasu  and  Matsumoto9  was  used  in  the 
MD/QCT  method  with  the  aforementioned  potential  energy 
surface. 

Energy  transfer  between  the  target  HC1  translational  and 
internal  energy  modes  was  investigated  and  it  was  found  that 
the  variation  of  the  inelastic  cross  section  has  a  negligible 
effect  on  the  transport  cross  section.  Hence,  MD/QCT  calcu¬ 
lations  are  able  to  predict  inelastic  cross  sections  and  allow 
for  a  more  accurate  estimation  of  the  viscosity  and  total  col¬ 
lision  cross  section  at  high  velocities.  The  latter  is  a  further 
test  of  the  adequacy  of  the  VHS  model  which  is  assumed  to 
be  correct  in  the  TCE  approach. 

Hypervelocity  collisions  such  as  O  +  HCl  occur  in  the 
jet-atmospheric  interaction  created  in  the  rarefied  space  en¬ 
vironment  when  a  reaction  control  system  (RCS)  is  posi¬ 
tioned  on  the  side  of  a  hypersonic  vehicle.  The  large  diffuse 
interaction  region  created  by  the  rarefied,  high  Mach  number 
chemically  reacting  flow  is  due  to  the  high  relative  velocities 
of  the  freestream  atomic  oxygen  and  the  HC1  plume  exhaust 
species.  The  accurate  modeling  of  the  interaction  region  is 
complicated  not  only  by  the  presence  of  hypervelocity  colli¬ 
sions,  but  also  by  the  complex  gas  dynamics.  The  evacuation 
of  a  thruster  with  concentration  approximately  four  orders  of 
magnitude  higher  than  the  low-pressure  space  environment 


into  the  diffuse  wind  creates  a  wide  range  of  characteristic 
length  scales  in  the  flow.  The  emphasis  of  earlier  work3'4  was 
to  establish  a  computational  approach  that  insured  an  accu¬ 
rate  DSMC  solution,  capturing  the  changing  flow  physics 
over  the  range  of  freestream  conditions  at  altitudes  of 
80-120  km  and  speeds  of  3-8  km/s.  The  bow-shock  inter¬ 
action  region  was  observed  to  have  continuum-like  features 
for  all  speeds  at  80  km,  and  the  strength  of  the  reaction  re¬ 
gion  in  terms  of  the  number  of  chemical  reactions  between 
freestream  and  thruster  species  was  found  to  scale  with  ve¬ 
hicle  velocity,  as  would  be  expected. 

In  this  work  we  propose  an  approach  for  modeling  hy¬ 
pervelocity  chemical  reactions  between  thruster  and  atmo¬ 
spheric  species  in  a  manner  suitable  for  DSMC  simulations. 
Since  HC1  is  present  in  high  concentration  in  the  thruster 
plume,  and  chemically  reactive  atomic  oxygen  is  available  as 
a  large  percentage  of  the  freestream,  accurate  modeling  of 
the  O  +  HCl  reaction  at  hypervelocity  conditions  is  crucial  to 
predicting  the  occurrence  of  OH,  an  important  radiator.  In 
fact,  at  altitudes  of  120  km  or  higher,  O  exchange  with  HC1 

was  found  to  contribute  more  than  70%  of  the  OH 

-2 

produced.  Although  other  plume  fragments  will  also  exist  in 
the  flow,  the  comparison  between  TCE  and  ab  initio 
MD/QCT  calculations  for  O  +  HCl  will  provide  an  assess¬ 
ment  of  the  overall  feasibility  of  using  TCE  to  model  similar 
hypervelocity  collisions. 

The  outline  of  the  remainder  of  this  paper  is  as  follows: 
Sec.  II  describes  some  problems  with  using  the  TCE  model 
for  hypervelocity  collisions  in  the  jet-atmospheric  interaction 
and  the  need  for  accurate  high  velocity  collision  cross  sec¬ 
tions.  The  MD/QCT  modeling  of  the  reaction  and  collision 
cross  sections  is  presented  in  Sec.  III.  In  Sec.  IV,  we  discuss 
the  specific  DSMC  numerical  parameters  used  for  the  single 
freestream  condition  examined  here:  120  km,  5  km/s.  Fi¬ 
nally,  in  Sec.  V,  we  present  detailed  comparison  of  the 
MD/QCT  results  with  the  TCE  model  as  well  as  more  accu¬ 
rate  high  temperature  O  +  HCl  viscosity  values  and  compare 
the  changes  observed  in  the  DSMC  simulation  for  the  differ¬ 
ent  chemical  reaction  parameters  for  the  O+HCl  reaction. 
Although  there  are  multiple  chemical  reactions  that  contrib¬ 
ute  to  the  production  of  OH,  we  limit  our  discussion  to  the 
O  +  HCl  reaction,  since  that  is  the  most  energetically  favor¬ 
able  reaction  path. 

II.  ISSUES  WITH  USING  THE  TCE  MODEL 

In  our  previous  work,  the  Arrhenius  parameters  from 
Mahmud  et  al.,1  A  =  5.6X  10-27  m3/ molecule/ s,  n=2.87,  and 
£,,=2.44  X  1  (A20  J,  were  used  for  TCE  calculations  for  the 
0+  HCI  — '  OH  +  CI  reaction,  but  in  this  work,  Arrhenius  pa¬ 
rameters  based  on  the  recent  calculations  of  Xie  et  al.10 
(shown  in  Table  I)  were  used  and  compared.  The  entire  set  of 
chemical  reactions  used  in  the  DSMC  simulations  can  be 
found  in  earlier  work.  None  of  the  other  reactions  were 
found  to  be  important  to  the  overall  flow  field  characteristics 
or  to  the  production  of  OH,  so  the  TCE  model  was  used  to 
calculate  all  other  reactions.  A  schematic  of  the  flow  geom¬ 
etry  for  the  RCS  plume-atmospheric  interaction  is  given  in 
Fig.  1.  We  will  consider  a  generic  RCS-vehicle  geometry 
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TABLE  I.  Freestream-plume  species  reactions  for  OH  production  reactions 
used  in  the  TCE  model. 


Reaction 

A,  m3/s 

n 

E„,  X  10“19  J 

h2o+n2 

^oh+h+n2 

5.81  X10"15 

0.00 

7.314 

H20  +  02 

->oh+h+o2 

1.13  X  10“7 

-1.31 

8.197 

h2o+o 

->OH+H+0 

1.13  X  10“7 

-1.31 

8.197 

h2o+o 

->OH+OH 

1.13  X  10“16 

0.00 

1.275 

h+o2^ 

OH+O 

1.66  X  10“16 

0.00 

1.061 

o+h2^ 

OH+H 

3.12X  10-16 

0.00 

0.952 

OH+C1- 

^O+HCl 

3.10X  10-27 

2.91 

0.070 

O  +  HCl- 

->OH+Cl  (Mahmud)a 

5.60  X  10-27 

2.87 

0.244 

O  +  HCl- 

^OH+Cl  (Xie)b 

1.70  X  10-22 

1.485 

0.408 

“Reference  7. 
bReference  10. 


with  freestream  conditions  corresponding  to  those  of  high 
altitude,  hypersonic  flight  for  a  lateral  side  jet  thrusting  per¬ 
pendicular  to  the  rocket  velocity  vector. 

To  illustrate  the  issues  involved  in  the  previous  calcula¬ 
tions,  Fig.  2  shows  the  spatial  distribution  of  the  TCE  reac¬ 
tion  probability  for  O  +  HCl— ►  OH + Cl  in  the  X-Y  plane  of  a 
sample  calculation  using  the  extrapolated  Mahmud  rates.7 
For  an  altitude  of  120  km  and  a  freestream  velocity  of 
5  km/s,  the  typical  interaction  region  ahead  of  the  thruster 
jet  is  seen.  It  is  observed  that  the  TCE  reaction  probabilities 
for  this  reaction  are  greater  than  unity  in  the  interaction  re¬ 
gion.  This  behavior  is  due  to  the  extrapolations  to  much 
higher  temperature/energy  regime  than  that  for  which  the 
Mahmud  parameters  are  valid,  as  was  mentioned  above.  If 
one  could  measure  reaction  rates  for  much  higher  tempera¬ 
tures,  it  is  extremely  unlikely  that  the  rate  coefficient  would 
continue  to  increase  with  the  same  steep  temperature  depen¬ 
dence.  Figure  3  shows  the  distribution  of  TCE  reaction  prob¬ 
ability  for  all  the  0+  HCI  — >OH  +  CI  reactions  in  the  entire 
domain,  showing  that  more  than  10%  of  the  collisions  have  a 
reaction  probability  greater  than  unity.  Since  a  probability 
greater  than  one  has  no  meaning  in  a  statistical  calculation, 
the  specific  implementation  used  in  Ref.  3  was  to  artificially 
set  these  probabilities  equal  to  unity.  This  limit  causes  a 
smaller  number  of  reactions  to  occur  in  the  simulation  com¬ 
pared  to  that  predicted  by  the  Arrhenius  rate. 

Another  question  that  affects  the  reaction  probabilities 
using  either  the  TCE  or  MD/QCT  model  is  the  assumed 
value  of  the  collision  cross  section.  The  definition  of  a  local 


x 


FIG.  2.  Reaction  probability  distribution  in  the  xy  plane  for  O  +  HCl 
— >OH+Cl  using  the  TCE  model  (rate  constant  of  Mahmud  et  al.,  Ref.  7)  for 
an  altitude  of  120  km  and  a  freestream  velocity  5  km/s. 


collision  frequency  based  on  cell  number  density  and  the 
velocity-dependent  collision  cross  section  for  any  sampled 
pair  of  molecules  is  fundamental  to  the  DSMC  method.  The 
collision  cross  section  is  most  commonly  defined  by  the  vari¬ 
able  hard  sphere  (VHS)  model,6  where  the  reference  diam¬ 
eter  and  velocity-dependence  are  defined  by  parameters  that 
reproduce  a  given  temperature-dependence  of  the  viscosity. 
There  are  large  uncertainties  in  the  values  of  the  VHS  model 
parameters,  which  are  obtained  from  relatively  low  tempera¬ 
ture  flows,  for  high  energies.  An  improved  approach  for  ob¬ 
taining  appropriate  high-energy  VHS  parameters  is  presented 
in  the  following  section. 
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FIG.  1.  Schematic  of  the  flow. 
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FIG.  3.  Reaction  probability  distribution  for  O+HCl— >  OH + Cl  using  the 
TCE  model  (rate  constant  of  Mahmud  et  al.,  Ref.  7)  for  an  altitude  of 
120  km  and  a  freestream  velocity  5  km/s. 
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Since  the  TCE  model  reproduces  the  prescribed  reaction 
rates  (under  equilibrium  conditions)  by  definition,  a  different 
value  in  the  VHS  collision  cross  section  does  not  affect  the 
reaction  rate  produced  in  a  computational  cell.  A  lower  value 
of  VHS  collision  cross  section,  for  instance,  would  increase 
the  reaction  probability,  but  then  would  correspondingly  de¬ 
crease  the  total  number  of  collisions  in  the  cell,  thus  produc¬ 
ing  the  same  overall  rate  of  reactions.  Hence,  the  assumed 
value  for  the  VHS  collision  cross  section  will  not  change  the 
number  of  reactions  in  a  simulation  using  the  TCE  model. 
However,  the  collision  cross  section  values  are  very  impor¬ 
tant  in  defining  the  reaction  probabilities  resulting  from  the 
MD/QCT  calculations.  These  probabilities  are  used  directly 
as  a  lookup  table  during  the  DSMC  simulation.  In  addition, 
comparisons  between  computed  reaction  probabilities  result¬ 
ing  from  MD/QCT  and  resulting  from  TCE  will  be  affected 
by  the  collision  cross  section  used. 


III.  MD/QCT  MODELING  OF  THE  REACTION 
PROBABILITY  AND  COLLISION  CROSS  SECTION 


The  MD/QCT  reaction  probability  is  the  ratio  of  the  re¬ 
action  cross  section  obtained  by  the  QCT-IEQMT  method  to 
the  total  collision  cross  section.  An  efficient  Monte  Carlo 
numerical  procedure  is  utilized  to  evaluate  the  multidimen¬ 
sional  integrals  associated  with  averaging  of  the  scattering 
geometric  properties  before  a  collision.  The  Monte  Carlo 
method  converges  at  a  rate  that  is  independent  of  the  dimen¬ 
sionality  of  the  integral.  In  order  to  utilize  the  Monte  Carlo 
method,  the  specific  initial  states  of  each  of  the  trajectories 
that  correspond  to  the  desired  reaction  cross  section,  such  as 
< rr(v,J,Ea )  or  (Tr(Eml,Elr),  where  v  is  the  vibrational  quan¬ 
tum  number,  J  is  the  rotational  quantum  number,  and  Ea  is 
the  translational  energy,  are  needed.  Microcanonical  sam¬ 
pling  enables  one  to  efficiently  sample  the  initial  states  of  the 
target  molecule  of  equal  energy  to  obtain  the  trajectory  initial 

13 

conditions  necessary  to  perform  the  MD/QCT  calculations. 
For  the  DSMC  simulations,  the  reaction  probability  of  a 
chemical  reaction  as  a  function  of  translational  energy,  Ea, 
and  molecular  internal  energy,  E)nt,  is  used.  The  reaction 
cross  section,  cr,.,  is  obtained  by  evaluating  whether  a  O 
+  HC1  collision  ends  in  a  chemical  reaction  by 


A.  MD/QCT  reaction  probability 

2 

A  parallel  MD/QCT  code"  was  developed  to  calculate 
the  reaction  and  viscosity  cross  sections  using  the  QCT- 
intemal  energy  quantum  mechanical  threshold  (QCT- 
IEQMT)  method.11  For  the  C)+  HCI  — ' OH  +  CI,  the  ab  initio 
',4 "  potential  energy  surface  of  Ramachandran  and  Peterson 
(RP)  (Ref.  8)  was  utilized.  The  RP  PES  is  accurate  for  total 
energies  less  than  40  kcal/mol.  However,  for  energies  higher 
than  40  kcal/mol,  the  reaction  of  O  +  HCl  may  have  two  pos¬ 
sible  outcomes. 


0(3P)  +  HC1(  ‘S+)  ->  OH(2II)  +  C1(2P)  (2) 


^C10(2ri)  +  H(2S).  (3) 


The  second  channel  forms  the  CIO  system,  a  stable  radical, 
but  a  potential  for  this  channel  is  not  available.  Xie  et  al.10 
estimated  the  contribution  from  this  channel  and  found  that 
at  3200  K  (a  maximum  total  energy  of  about  4.42  X  10-19  J), 
the  rate  constant  for  this  channel  was  roughly  two  orders  of 
magnitude  smaller  than  the  total  rate  constant.  At  120  km 
altitude  with  a  freestream  velocity  of  5  km/s,  most  of  the 
O  +  HCl  collisions  have  the  energy  less  than  5  X  10-19  I. 
Thus,  in  this  work,  we  neglect  this  channel.  The  RP  surface 
was  also  utilized  by  Xie  et  al.,  who  computed  both  quantum 
and  QCT  reaction  cross  sections  for  a  range  of  energies. I- 
However,  it  was  necessary  for  us  to  perform  our  own 
MD/QCT  calculations  since  a  large  set  of  probabilities  are 
required  for  the  DSMC  flow  field  simulation,  and  we  require 
cross  sections  for  high  total  collision  energy  values.  In  addi¬ 
tion,  the  DSMC  model  used  here  assumes  that  the  vibrational 
and  rotational  energies  are  continuous,  thus  requiring  a  table 
of  probabilities  for  a  range  of  internal  energies,  Em,  while 
Xie  et  al.  computed  reaction  cross  sections  for  the  vibra¬ 
tional  quantum  number,  v  =  0,  1 ,  and  2  states  only.  It  will  be 
shown  in  Sec.  IV  that  the  MD/QCT  results  presented  here 

agree  well  with  those  of  Xie  et  al.  for  the  energy  range 
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considered  in  that  work. 


=  7 Tb 


2 

max 


Nr 

Nt' 


(4) 


where  b  is  the  impact  parameter,  /;max  is  the  maximum  im¬ 
pact  parameter  beyond  which  a  chemical  reaction  does  not 
occur,  Nr  is  the  number  of  trajectories  that  result  in  a  reac¬ 
tion,  and  NT  is  the  number  of  total  trajectories.  In  previous 
work,-  the  reaction  probability,  Pn  for  another  chemical  re¬ 
action  was  assumed  to  be 


1  r  w/ 

°VHS 

For  the  O  +  HCl  case,  however,  it  will  be  shown  that  the 
MD/QCT  collision  cross  section  is  different  than  that  given 
by  the  VHS  model  and  will  be  discussed  further  in  the  next 
subsection. 


B.  Collision  cross  section 


As  mentioned  above,  in  our  previous  work,  the  variable 
hard  sphere  (VHS)  model  with  Bird  values  was  used  for  the 
collision  cross  section.  The  VHS  cross  section  is  defined  by 
crVHS=Trc/2  with6 


(2  kTieir 

(mrg2Y°T(2  -  co) 


1/2 


(6) 


where  d  is  the  diameter  of  a  molecule,  mr  is  the  molecular 
reduced  mass,  g  is  the  relative  velocity,  and  oj  is  the  viscosity 
index  (a>=v-0.5,  v  is  the  temperature  exponent  of  the  vis¬ 
cosity  coefficient).  «=0.375,  t/ref=4.38  A,  and  2j.ef=273  K 
for  the  O  +  HCl  reaction  in  Ref.  5.  The  reference  diameters, 
dre f,  are  calculated  from  the  viscosity  data, 

(  MmkTJir)m  \1/2 

ref  \  4(4  -  2<u)(6  -  2w)yu,ref/  ’  1  J 


where  /jiKi  is  the  reference  viscosity  coefficient,  and  m  is  the 
mass  of  the  molecule. 
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FIG.  4.  Comparison  of  distributions  of  the  deflection  angle  \  f°r  O+HCl 
collisions  predicted  by  the  VHS,  VSS,  LDF,  and  MD/QCT.  The  MD/QCT 
results  were  calculated  for  £int=2.0  X  10-19  J  and  a  relative  velocity  of 
5  km/s. 


One  could  also  attempt  to  calculate  the  total  collision 
cross  section  for  the  RP  potential  surface,  but  the  unbounded 
nature  of  a  classical  total  cross  section  for  a  realistic  potential 
comes  into  play.  To  illustrate,  one  may  use  the  MD/QCT 
method  to  compute 

°VHS  =  ^max^-  >  (8) 

where  Nc  is  the  number  of  trajectories  that  result  in  a  colli¬ 
sion  and  Nt  is  the  number  of  total  trajectories.  However,  in 
order  to  decide  whether  a  collision  has  occurred,  the  cutoff 
deflection  angle,  x ,  is  needed,  and  ctVhs  was  found  to  be 
strongly  dependent  on  the  cutoff  angle  x  f°r  the  O  +  HCl 
system. 

Figure  4  shows  a  comparison  of  the  distribution  of  the 
deflection  angle  x  f°r  O+HCl  collisions  predicted  by  the 
VHS,  and  variable  soft  sphere  (VSS)  (Ref.  14)  models,  the 
linear  deflection  function  (LDF),  and  MD/QCT  results.  The 
MD/QCT  results  were  calculated  at  £’int=2.0  X  10-19  J  and 
the  relative  velocity  of  5  km/s.  This  HC1  internal  energy  is 
typical  for  the  flow  field  freestream  conditions  at  120  km  and 
velocity  of  5  km/s.  For  the  VHS  and  VSS  models,  the  de¬ 
flection  functions  are  given  by 


XLDF(b)=m{l-b/dLDF),  (11) 

where  a  value  of  dLDF= 2.63  A  was  chosen  to  match  with  the 
MD/QCT  results.  It  is  shown  that  the  LDF  model  agrees 
better  with  the  MD/QCT  results  than  the  VHS  and  VSS  mod¬ 
els.  However,  the  discrepancy  increases  between  the 
MD/QCT  results  and  all  of  the  models  for  smaller  deflection 
angles.  The  question  of  whether  the  DSMC  collision  cross 
sections  should  be  parametrized  by  the  collider  pair  internal 
energy  has  been  discussed,15  but  will  not  be  considered  here. 
For  higher  internal  and  translational  energies,  the  cutoff 
angle  is  more  important,  and  it  is  difficult  to  define  as  may  be 
seen  in  the  figure.  Therefore,  the  approach  that  Tokumasu9 
used  is  utilized  to  calculate  the  viscosity  cross  section. 

For  implementation  in  DSMC,  it  is  more  consistent  to 
obtain  a  transport-based  collision  cross  section  (usually  vis¬ 
cosity)  from  which  to  convert  a  reaction  cross  section  to  a 
reaction  probability,  since  the  DSMC  method  already  uses 
the  VHS  model  viscosity-based  cross  section  as  the  basis  for 
computing  the  number  of  collisions  per  cell  per  time  step  and 
as  the  basis  for  collision  dynamics  (hard-sphere  isotropic 
scattering).  In  this  spirit,  the  work  of  Tokumasu  and 
Matsumoto9  demonstrated  the  calculation  of  an  accurate  vis¬ 
cosity  cross  section  for  a  given  potential  using  a  Monte  Carlo 
method  for  integration.  While  the  integration  over  the  impact 
parameter  becomes  infinite  for  the  total  collision  cross  sec¬ 
tion,  those  for  classical  momentum  and  energy  transfer  cross 
sections  are  finite.  Therefore,  the  viscosity  cross  section  was 
calculated  first,  and  converted  to  the  equivalent  VHS  colli¬ 
sion  cross  section  for  each  collision  velocity.9 

The  viscosity  cross  section,  cr^,  is  calculated  by  the 
Monte  Carlo  evaluation  of  an  integral  given  by  Ref.  16 

<V,md  =  J  (  ^  sin2  x  +  ^(A<?int)2  -  ~(Aeint)2  sin2  X^dr, 

(12) 

where  g  is  the  dimensionless  relative  velocity  that  is  defined 
as  \mr/kTg,  eim  is  the  dimensionless  internal  energy  that  is 
nondimensionalized  by  kT,  and  A<?int  is  the  change  of  the 
dimensionless  HC1  internal  energy  before  and  after  a  colli¬ 
sion.  The  quantities  g,  and  Aeint  are  obtained  from  the 
MD/QCT  calculations  similar  to  those  performed  for  the  re¬ 
action  probability  except  that  the  sampling  of  initial  trajec¬ 
tory  conditions  corresponds  to  the  geometric  precollisional 
conditions  specified  by  Aeint.  The  MD/QCT  viscosity  cross 
section  cr^  MD  converges  when  the  maximum  impact  param¬ 
eter  is  selected  so  that  the  effect  of  the  interaction  potential  is 
negligible.  The  integral  f()dr  specifies 


*vhs  =  2cos  \bld). 

(9) 

A'vss  =  2cos_1(/?/g01/“, 

(10) 

where  a  is  an  adjustable  parameter  between  1  and  2,  and  a 
value  of  1.59  was  used  in  the  figure.5  For  this  comparison,  a 
diameter,  d,  of  2.35  A  was  used.  The  LDF  model  gives  the 
deflection  angle  as 


f  fb  max  f"  f2’r  C* 

\°dr=\  I  I  I 
J  J  o  Jo  Jo  J  c 


0  JO  Jo  JO 


dip 

77 


d(p 

277 


1 

sin  Xdx 


X[2irbdb\ , 


03) 


where  ift  is  the  initial  orientation  of  molecules  and  <f>  is  the 
initial  direction  of  the  rotational  vector.  In  addition,  the  vis- 

17 

cosity  coefficient,  //,  is  obtained  by 
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TABLE  II.  Freestream  parameters. 


Parameter 

120  km 

Temperature,  K 

354 

Number  density,  molecule /m3 

4.73  X1017 

02  mole  fraction,  % 

9 

N2  mole  fraction,  % 

73 

O  mole  fraction,  % 

18 

]_ 


8  1 
5  'j2TrmrkT  6int 


dJ 


xj  dE^J^(2J+  l)exp' 


kT 


(14) 


where  Qmt  is  the  partition  function  for  internal  modes.  The 
equivalent  VHS  collision  cross  section  is  crVHS  MD=  7rc/yHS, 
and  the  diameter  for  the  collision  cross  section  can  be  ob¬ 
tained  by 

^VHS-  \  „4  ■  (15) 

V  7 Tg 


IV.  DSMC  MODELS,  NUMERICAL  PARAMETERS, 

AND  FREESTREAM  CONDITIONS 

Figure  1  shows  the  flow  geometry  used  in  the  modeling 
of  the  interaction  of  the  atmosphere  and  the  RCS  jet.  A  small 
rocket  is  modeled  as  a  blunted  cone  cylinder  and  a  thruster 
positioned  on  the  cylinder  right  after  the  cone-cylinder  junc¬ 
tion.  The  radius  of  the  cylinder  is  0.2  m,  the  length  from  the 
head  of  the  cone  to  the  nozzle  exit  is  2  m,  and  the  angle  of 
attack  is  zero.  The  freestream  parameters  at  120  km  altitude 
are  listed  in  Table  II.  The  species  mole  fractions  at  the  nozzle 
exit  is  also  listed  in  Table  III.  A  starting  surface  was  obtained 
from  the  axisymmetric  plume  core-flow  DSMC  simulations. 
The  plume  core-flow  DSMC  simulations  were  performed 
with  the  nonuniform  nozzle  exit  condition.  '4  The  density 
isolines  of  about  6X  1021  molecules/m3  were  taken  for  the 
starting  surface  of  a  60  Ibf  (270  N)  thruster.  The  starting 
surface  is  an  oval  shape  with  approximate  x  and  y  dimen¬ 
sions  of  0.3  and  0.5  m. 

The  entire  set  of  the  chemical  reactions  between  the 
thruster  side  jet  and  plume-atmospheric  species  includes  re¬ 
actions  between  oxygen  and  nitrogen  species  as  well  as  ones 
that  either  produce  or  consume  OH  (see  Table  I).  The  first 


TABLE  III.  Nozzle  exit  conditions. 


Species 

Mole  fraction 

h2o 

25% 

n 

o 

5% 

CO 

23% 

HC1 

14% 

n2 

14% 

h2 

19% 

FIG.  5.  Reaction  rate  constant  for  O+HCl— »OH+Cl  as  a  function  of  tem¬ 
perature.  The  MD/QCT(RP)  is  calculated  in  this  work. 


three  reactions  in  Table  I  involve  the  dissociation  of  water  by 
freestream  constituents,  N2,  02,  and  O.  Note  that  the  relative 
importance  of  these  three  reactions  will  change  with  altitude. 
The  last  four  reactions  are  exchange  reactions  between  the 
freestream  and  plume  constituents.  The  following  two  ex¬ 
change  reactions,  0  +  H2  and  O  +  HCl,  are  potentially  impor¬ 
tant  because  both  H2  and  HC1  are  present  in  relatively  large 
mole  fractions  of  divert  solid  propellant  motors. 

The  three-dimensional  DSMC  calculations  were  imple¬ 
mented  in  the  Statistical  Modeling  In  Low-density  Environ¬ 
ment  (SMILE)  (Ref.  18)  computational  tool.  The  gas  is  con¬ 
sidered  a  14-species  reacting  mixture.  Approximately  18 
million  molecules  were  simulated  in  the  computational  do¬ 
main.  Separate  grids  were  used  for  collisions  and  macropa¬ 
rameters  adaptive  to  flow  gradients,  and  the  total  number  of 
collision  cells  was  approximately  5  million.  The  total  num¬ 
ber  of  time  steps  was  about  100  000  with  a  time  step  of 
2.0  X  10-7  s,  and  macroparameter  sampling  was  started  after 
20  000  time  steps,  a  time  sufficient  to  reach  the  steady  state. 
These  parameters  were  chosen  to  minimize  the  statistical  de¬ 
pendence  between  simulated  particles,  remove  the  grid  de¬ 
pendence  of  the  results,  and  furnish  sufficient  spatial  resolu¬ 
tion  of  the  boundary  layer  along  the  rocket. 

The  major  frequency  scheme18  is  employed  for  modeling 
the  molecular  collision  frequency  to  enhance  the  statistical 
representation.  For  trace  species,  a  weighting  factor  of  ap¬ 
proximately  0.01  was  used.  The  VHS  model  was  used  for 
modeling  nonreactive  interactions  between  particles.  The 
Borgnakke -Larsen1 7  (BL)  model  (see  also  Ref.  20)  with 
temperature-dependent  rotational  and  vibrational  relaxation 
numbers  was  chosen  for  modeling  rotation-translation  and 
vibration-translation  energy  transfer.  The  effective  degree  of 
freedom  for  vibration  was  modeled  to  be  temperature- 
dependent.  The  gas-surface  interaction  was  modeled  using  a 
diffuse  model  with  total  energy  and  momentum  accommoda¬ 
tion  with  a  surface  wall  temperature  of  300  K. 
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V.  RESULTS  AND  DISCUSSION 
A.  Rate  constant  for  0  +  HCI->OH  +  CI 


Although  the  DSMC  method  requires  the  probability  of 
a  reaction,  instead  of  the  rate  constant,  the  latter  quantity 
provides  a  consistency  check  for  the  MD/QCT  calculations 
presented  in  this  work.  Figure  5  presents  the  reaction  rate 
constant  calculated  using  a  parallel  MD/QCT  code  between 
1000  K  and  3000  K.  The  reaction  rate  may  be  calculated 
from  the  reaction  cross  section  by  averaging  the  target  mol¬ 
ecule  (HC1)  over  a  Maxwellian  equilibrium  distribution, 

V kTl 


Kf=f(T)Q~LQ 


irmrj 

f00  Cbm  ax 

xSH(2/+l)  Pr(v,J,Eu 

v  J  J 0  Jo 


b) 


X  exp(-  EtlJkT)2TrbdbEtrdEtr, 


(16) 


where  Qvib  and  Qmt  are  the  vibrational  and  rotational  parti¬ 
tion  function,  respectively,  and  for  O  +  HCl, 


f(T)  = 


_ 3 _ 

5  +  3  exp(-  228 IT)  +  exp(-  326 IT) 


(17) 


The  factor  f{T)  is  the  probability  that  the  target  molecular 
system  is  initially  on  one  of  the  three  electronic  surfaces  that 
allow  a  reaction  to  occur.  This  temperature  dependent  ex¬ 
pression  accounts  for  the  spin-orbit  splitting  of  the  overall 
triplet  reagents.  The  typical  number  of  trajectories  per  tem¬ 
perature  value  was  approximately  20  000.  Figure  5  shows 
that  the  rate  of  Mahmud  et  al.1  used  in  our  previous  DSMC 
calculations  is  higher  than  the  MD/QCT(RP)  rate  obtained 
with  the  RP  surface  in  this  work  and  the  improved  canonical 
variational  theory  (ICVT)  rate  constant  of  Xie  et  al.w  The 
MD/QCT(RP)  rate  is  slightly  lower  than  the  ICVT  rate  be¬ 
cause  only  the  iA"  surface  was  used  for  MD/QCT  calcula¬ 
tions,  and  the  3,4 '  surface  contribution  increases  for  higher 
temperature.  In  addition,  the  tunneling  effect,  not  included  in 
the  MD/QCT  calculations,  is  more  important  for  lower  tem¬ 
peratures. 


B.  Reaction  probability  for  O  +  HCl— >OH  +  CI 

Figure  6  shows  a  comparison  of  reaction  probabilities 
between  the  MD/QCT  and  TCE  (Ref.  6)  models  as  a  function 
of  reactant  internal  energy  at  3  and  5  km/s  relative  veloci¬ 
ties.  The  MD/QCT  reaction  probability  is  calculated  here  as 
the  reaction  cross  section  divided  by  the  VHS  cross  section 
using  Bird’s  values.  The  rate  of  Mahmud  et  al.  used  for  our 
previous  simulations  and  the  rate  of  Xie  et  al. 10  are  also 
shown  in  the  figure.  The  effective  number  of  degrees  of  free¬ 
dom  for  vibration  £„hci=  1  is  used,  and  the  VHS  model  with 
the  parameters  of  Ref.  5  of  <n= 0.375,  ofref=4.38  A,  and  Tief 
=  273  K  was  used  for  all  the  three  cases  in  this  subsection.  In 
the  next  subsection,  the  fidelity  of  these-specific  VHS  values 
used  will  be  discussed.  It  can  be  seen  that  the  TCE  probabili¬ 
ties  computed  from  the  rate  of  Mahmud  et  al.1  are  much 
higher  than  the  MD/QCT  or  TCE  model  using  the  Arrhenius 
parameters  derived  from  the  rate  coefficients  computed  by 
Xie  et  al. 10  for  both  relative  velocities.  This  result  is  consis- 


FIG.  6.  Comparison  between  MD/QCT  and  TCE  models:  Reaction  prob¬ 
abilities  for  O+HCl— » OH + Cl  as  a  function  of  the  reactant  internal  energy. 


tent  with  our  previous  calculations  in  which  the  probabilities 
for  the  O+HCl  reaction,  using  the  TCE  model  and  Mahmud 
extrapolated  rates,  were  found  to  be  unrealistically  high.  For 
the  3  km/s  case,  the  MD/QCT  model  predicts  slightly  higher 
reaction  probabilities  than  the  TCE  with  the  rate  calculations 
of  Xie  et  al.10  Nonetheless,  at  5  km/s,  the  discrepancy  be¬ 
comes  more  significant,  and  the  MD/QCT  model  predicts 
higher  probabilities  than  the  TCE  model. 

C.  Collision  cross  section  for  O  +  HCl 

In  the  previous  subsection,  the  VHS  collision  cross  sec¬ 
tion  with  <w=0.375,  t/ref=4.38  A  and  7j.ef=273  K  was  used. 
In  this  subsection,  the  preferred  values  for  VHS  for  O 
+  HC1  collisions  is  investigated.  In  our  earlier  reported 
results,3  4  the  viscosity  index  w  of  0.25  for  O  and  co  of  0.5  for 
HCl  were  used,  ’  i.e.,  the  HC1  molecule  was  treated  in  the 

Maxwell  model.  However,  although  the  viscosity  index  a>  of 
HCl  is  0.5  between  20  and  99  °C,  for  some  gases,  to  de¬ 
creases  as  temperature  increases."1  Therefore,  the  VHS  col¬ 
lision  cross  sections  obtained  from  a  Maxwell  molecule  may 
not  be  accurate  for  higher  temperatures. 

Let  us  consider  the  selection  of  two  parameters:  (1)  tem¬ 
perature  exponent  of  the  viscosity  coefficient,  v,  and  (2)  the 
reference  diameter,  dTef,  for  the  VHS  model.  The  parameters 

used  in  the  previous  DSMC  calculations  were  mostly  ob- 

21 

tained  from  viscosity  data  in  the  low  temperature  range." 
Figure  7  shows  the  variation  of  viscosity  coefficient  of  HCl 
from  100  to  5000  K.  In  the  figure,  the  data  of  Svehla"-  and 
the  data  (Bird)  in  Ref.  21  are  shown.  Svehla  estimated  the 
viscosities  using  the  Lennard-Jones  (12-6)  potentials.  From 
the  data  (Bird),  the  temperature  exponent  of  the  viscosity 
coefficient  of  HCl,  i'hci,  is  TO,  and  with  vHC1  of  1.0,  the 
viscosity  coefficient  is  higher  than  the  data  of  Svehla.  The 
temperature  exponent  of  the  viscosity  coefficient  of  HCl, 
Viicb  °f  0-65  agrees  well  with  the  data  of  Svehla  between 
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FIG.  7.  HC1  viscosity  coefficient,  /jl ,  for  high  temperatures.  Svehla’s  viscos¬ 
ity  coefficient  is  compared  with  the  Maxwell  model  (Bird,  j/=1.0).  vHCI  of 
0.65  gives  great  agreement  with  the  data  of  Svehla. 

1000  and  5000  K.  The  VHS  parameters  for  HC1  are  listed  in 
Table  IV.  It  is  seen  that  the  value  of  co  is  quite  different  for 
low  and  high  temperatures. 

Using  the  RP  3 A"  surface  in  the  MD/QCT  calculations, 
the  quantities  such  as  deflection  angle  and  the  change  in  HC1 
internal  energy  are  obtained  and  used  in  Eq.  (12)  to  calculate 
the  viscosity  cross  section.  The  change  in  HC1  internal  en¬ 
ergy,  eint,  is  obtained  by  analyzing  the  momentum  of  the 
postcollisional  HC1  species  in  each  trajectory.  To  account  for 
the  neglect  of  the  RP  3 A"  potential,  the  maximum  impact 
parameter  was  selected  such  that  the  rate  of  translational  en¬ 
ergy  transfer  was  less  than  10%,  and  the  MD/QCT  viscosity 
cross  section  cr^  MD  converged.  The  approach  that  Tokumasu 
used  is  this  work  to  obtain  both  the  O  +  HCl  viscosity  coef¬ 
ficients  as  well  as  the  VHS-equivalent  collision  cross  section. 

Figure  8  shows  a  comparison  of  the  viscosity  coeffi¬ 
cients  obtained  from  the  MD/QCT  results  and  the  O  and  HC1 
data  in  Ref.  21  (Bird)  and  data  of  Svehla.-”  To  compare  with 
the  MD/QCT  result  we  take  the  average  of  the  separate  O 
and  HC1  viscosity  coefficients.  The  MD/QCT  results  for  O 
+  HC1  are  believed  to  be  the  most  accurate  and  fall  above  or 
below  the  data  of  Refs.  21  (Bird)  and  22,  depending  on  the 
temperature.  The  agreement  between  the  MD/QCT  results 
and  the  viscosities  derived  from  the  low-temperature  data  of 
Ref.  21  is  good.  At  higher  temperatures,  the  difference  in  the 
MD/QCT  results  and  the  viscosities  derived  from  the  coeffi¬ 
cients  of  Svehla  is  attributed  to  the  simpler  potential  model 
used  in  the  latter  calculations. 

Figure  9  presents  a  comparison  of  the  collision  cross 


TABLE  IV.  VHS  parameters  for  HCI. 


dm{,  A 

Data 

CO 

Tef,  K 

Bird 

0.50 

5.76 

273 

Svehla 

0.15 

3.15 

5000 

FIG.  8.  O  and  HCI  mixture  viscosity  for  high  temperatures  calculated  by  the 
MD/QCT  method  [Eq.  (14)]. 


sections  between  the  MD/QCT  and  VHS  models  at  Emt 
=  0.5  X  10-19  and  2.0  X  10-19  J.  For  the  two  viscosity  param¬ 
eters  tested  in  the  VHS  model.  Fig.  9  shows  that  the  collision 
cross  sections  are  greater  with  high  temperature  data  of 
Svehla.  The  VHS  cross  section  with  wHC1=0. 15  and 
t/Hci,rL,f=5. 1 5  at  5000  K  (Svehla)  is  significantly  larger  than 
the  VHS  cross  section  with  wHC1=0.5  and  </hci, ref  =5-76  at 
273  K  (Bird)  by  more  than  50%.  The  MD/QCT  results  are 
between  the  VHS  cross  sections  obtained  from  the  data  of 
Ref.  5  (Bird)  and  those  of  Svehla.22  Therefore,  if  the 
MD/QCT  VHS  equivalent  collision  cross  section  in  Eq.  (15) 
is  used,  lower  reaction  probabilities  are  predicted  than  those 
that  would  be  obtained  from  the  VHS  cross  sections  using 
the  data  in  Ref.  5  (Bird).  Also,  it  was  found  that  the 
MD/QCT  VHS  equivalent  cross  section  does  not  change  sig- 


FIG.  9.  Comparison  of  collision  cross  sections  between  MD/QCT  and  VHS 
models  for  the  O  +  HCl  collision. 
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FIG.  10.  Comparison  of  inelastic  collision  cross  sections  for  the  O+HCl 
collision  at  £int=0.5X  10“19  and  Eint=2.0X  10"19  J. 


nificantly  with  changes  of  the  HC1  internal  energy.  Since  the 
MD/QCT  VHS  equivalent  cross  section  is  not  dependent  on 
the  HC1  internal  energy,  the  cross  sections  could  be  fit  to  the 
simpler  VHS  model.  The  parameters,  «=0.39  and 
t/ref.iooo  k=3-9  A  were  found  to  give  a  good  fit  of  the 
MD/QCT  cross  sections  to  the  VHS  form.  The  collision 
cross  sections  derived  from  these  parameters  will  be  desig¬ 
nated  as  crVHS  MD. 

The  inelastic  collision  cross  section,  crin,  is  defined  for 
collisions  that  do  not  react  but  produce  a  change  in  Elr  as 

of, r=  J  (A ea)2dr.  (18) 

The  quantity  Aetr  is  the  change  of  the  dimensionless  transla¬ 
tional  energy  before  and  after  a  collision  and  was  calculated 
from  the  MD/QCT  pre  and  post  trajectory  HC1  conjugate 
momenta  for  nonreactive  trajectories  as 

g'2  g2 

Ae,r=  '—  +  eCM~  ~  (19) 

where  eCM  is  the  dimensionless  center-of-mass  translational 
energy,  and  the  prime  indicates  postcollisional  quantities. 
Figure  10  shows  the  inelastic  collision  cross  section  as  a 
function  of  the  relative  velocity.  As  shown  in  the  figure,  the 
inelastic  cross  section  decreases  as  the  translational  energy 
increases.  It  is  also  found  that  the  inelastic  cross  section  is  a 
function  of  both  the  translational  and  internal  energy.  How¬ 
ever,  because  the  change  of  crin  is  insignificant  compared  to 
the  VHS -equivalent  collision  cross  section,  crVHS  MD  does  not 
show  a  strong  dependence  on  internal  energy. 

Figures  11  and  12  show  a  comparison  of  the  reaction 
probability  for  0+  HCI  —  OH  +  CI  calculated  by  using  the  re¬ 
action  cross  section  divided  by  either  the  VHS  model  (Bird) 
or  crVHSMD  at  3  km/s  and  5  krn/s,  respectively.  The  VHS 
model  with  co=039  and  drefjooo  k=3-9  A  was  used  to  repre¬ 
sent  crVHS  MD  because  the  MD/QCT  collision  cross  section 


FIG.  11.  Comparison  of  the  reaction  probability  calculated  by  using  the 
reaction  cross  section  divided  by  either  the  VHS  model  (Bird)  or  the 
MD/QCT  collision  cross  section  (trVHS  MD):  Reaction  probabilities  as  a  func¬ 
tion  of  the  reactant  internal  energy  at  3  km/s. 

results  agree  well  with  the  VHS  model  with  those  parameters 
as  shown  in  Fig.  9.  Both  the  MD/QCT  and  TCE  reaction 
probabilities  using  crVHSMD  are  lower  than  those  with  the 
VHS  cross  section  (Bird,  toHC1=0.5)  because  crVHS  MD  are 
greater  than  the  VHS  cross  sections  (Bird,  «HC]  =  0.5).  How¬ 
ever,  the  change  in  the  TCE  model  is  smaller  than  that  in  the 
MD/QCT  model. 

Figure  13  shows  computed  MD/QCT  reaction  probabil¬ 
ity  values  (using  ovhs.md)  as  a  function  of  total  collision 
energies.  The  filled  and  open  symbols  show  the  reaction 
probabilities  for  a  given  collision  velocity  with  varying  in¬ 
ternal  energy  and  with  a  fixed  internal  energy  with  varying 
collision  velocity,  respectively.  The  results  show  no  signifi¬ 
cant  favoring  effect  of  internal  energy  on  the  reaction  prob- 


HCI  internal  energy,  J 

FIG.  12.  Comparison  of  the  reaction  probability  calculated  by  using  the 
reaction  cross  section  divided  by  either  the  VHS  model  (Bird)  or  the 
MD/QCT  collision  cross  section  (trVHS  MD):  Reaction  probabilities  as  a  func¬ 
tion  of  the  reactant  internal  energy  at  5  km/s. 
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FIG.  13.  Reaction  probability  as  a  function  of  total  collision  energy  for 
different  collision  velocities  and  HC1  internal  energies.  Note  that  3,  5, 
7  km/s  corresponds  to  8.3  X  10-2°,  2.3  X  10-19,  4.5  X  10-19  J.  respectively. 


•  23 

ability.  This  result  is  somewhat  surprising  in  that  Aoiz  et  al. 
found  significant  vibrational  favoring  effects  in  their  QCT 
results  for  this  reaction.  However  a  less  accurate  PES  was 
used  in  that  study. 

Figure  14  presents  the  distributions  of  the  Aeint  at 
5  km/s  for  O  +  HCl  collision  for  two  initial  internal  energies. 
Aeint  was  calculated  by  Aeint=e'nt-eint.  Figure  14  shows  that 
as  the  initial  HC1  internal  energy  increases,  the  Aeinl  distri¬ 
bution  spreads.  In  other  words,  the  inelastic  cross  section 
increases,  as  the  initial  internal  energy  increases,  as  is  gen¬ 
erally  expected  based  on  simple  models  (see,  e.g..  Ref.  24). 
At  £int=0.5  X  1CT19  J,  energy  is  mostly  transferred  from  the 
translational  to  HC1  internal  modes.  Figure  15  shows  the 
distributions  of  the  Aeint  at  £'int=2.0  X  1CT19  J  for  different 


1C)'3  ■ 

-1.5  -1  -0.5 


FIG.  15.  Distribution  of  the  Aeint  for  £int=2.0X10  19  J  at  3  km/s  (top), 
5  km/s  (middle),  and  8  km/s  (bottom)  for  the  O+HCl  collision. 


collisional  translational  energies.  At  3  km/s,  the  energy  is 
transferred  from  the  HC1  internal  to  the  translational  mode. 
However,  as  the  initial  translational  energy  increases,  the  en¬ 
ergy  shifts  from  the  translational  to  internal  modes.  In  gen¬ 
eral  the  results  show  a  tendency  toward  energy  equilibration 
between  modes.  When  Etr >Eint,  the  Eint  gains  and  when 
£’int>£’tn  the  Etr  gains. 


FIG.  14.  Distribution  of  the  Aeint  for  Emt=0.5  X  10  19  J  (top)  and  £’int=2.0 
X  10“ 19  J  (bottom)  at  5  km/s  for  the  O+HCl  collision. 


FIG.  16.  N2  number  density  (molecule/ m3)  contours  with  the  rate  of  Xie  et 
al.  (TCE)  at  120  km  altitude  for  freestream  velocity  of  5  km/s.  Area  shown 
is  7.4  X  7.5  m. 
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FIG.  17.  Translational  temperature  contours  with  the  rate  of  Xie  et  al. 
(TCE)  at  120  km  altitude  for  freestream  velocity  of  5  km/s.  Area  shown  is 
7.4  X  7.5  m. 
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FIG.  18.  Mach  number  contours  with  the  rate  of  Xie  et  al.  (TCE)  at  120  km 
altitude  for  freestream  velocity  of  5  km/s.  Area  shown  is  7.4  X  7.5  m. 


D.  DSMC  calculations 

First,  the  TCE  model  with  the  rate  constant  of  Xie  et 
al.w  for  the  C)  +  HCI  — >OH  +  CI  reaction  was  utilized  for  the 
DSMC  calculation  at  120  km  altitude  for  freestream  velocity 
of  5  km/s.  For  the  TCE  model,  the  VHS  collision  cross  sec¬ 
tion  for  O  +  HCl  with  o>=0.375,  c/ref=4.38  A  and  7j.ef 
=  273  K  was  used.  Figure  16  shows  the  N2  number  density 
(molecule/m3)  contour  for  this  case.  The  maximum  N2  num¬ 
ber  density  is  about  1.2  X  1018  rrT3  located  at  1-2  m  ahead  of 
the  rocket  cone.  It  was  found  that  the  change  of  chemical 
reaction  models  and  O  +  HCl  collision  cross  section  do  not 
affect  the  overall  flow  field  because  collisions  between  O  and 
HCI  are  not  the  main  collision  pairs  (O  is  18%  of  the 
freestream  and  HCI  is  14%  of  the  side  jet).  Figure  17  shows 
the  overall  translational  temperature  contour  with  the  rate 
constant  of  Xie  et  al.  for  the  O  +  HCl— >  OH + Cl  reaction  at 
120  km  altitude  for  freestream  velocity  of  5  km/s.  Similar  to 
the  total  number  density,  the  overall  translational  temperature 
profile  does  not  change  for  the  different  O  +  HCl  reaction 
probabilities.  While  the  high  number  density  shock  region  is 
located  between  x=-2  and  -1  m,  the  high  temperature  shock 
region  is  shifted  more  forward  than  the  high  number  density 
shock  region.  Figure  18  shows  the  Mach  number  contours  at 
120  km  altitude  for  a  freestream  velocity  of  5  km/s.  The 
high  temperature  shock  region  coincides  with  the  low  Mach 
number  region. 

Second,  the  DSMC  calculations  were  implemented  for 
the  three  O  +  HCl  chemical  reaction  models  at  120  km  alti¬ 
tude  for  freestream  velocity  of  5  km/s.  In  the  discussion 
below,  we  designate  the  TCE  model  used  with  the  previous 
rate  constant7  as  case  (1);  the  TCE  model  using  the  rate 
constant  of  Xie  et  al.10  as  case  (2);  and,  the  MD/QCT  reac¬ 
tion  probability  as  case  (3).  For  case  (3),  the  MD/QCT  reac¬ 
tion  probability  used  in  this  subsection  is  the  ratio  of  the 
MD/QCT  reaction  cross  section  to  the  <tvhs  MD,  and  the  com¬ 
puted  MD/QCT  reaction  probabilities  were  tabulated  and 
used  in  the  DSMC  simulations.  For  cases  (1)  and  (2),  the 


O+HCl  collision  cross  section  used  in  calculating  the  elastic 
cross  section  for  pair-selection  and  elastic  collision  outcomes 
was  the  VHS  model  with  the  original  parameters  of  u> 
=  0.375  and  <7ref,273  k=4-38  A.  For  case  (3),  the  VHS  model 
with  the  modified  parameters,  <u=0.39  and  dle f  1000  k 
=  3.9  A  was  used.  The  change  in  collision  cross  sections 
changes  the  number  of  collisions  between  O  and  HCI  mol¬ 
ecules  which  in  turn,  affects  the  OH  production  rate. 

Figure  19  shows  the  OH  number  density  contours  with  a 
maximum  OH  number  density  of  about  2  X  1016  rrT3  with 
the  case  (2)  chemical  reaction  rate.  In  Figs.  20  and  21,  a 
quantitative  comparison  between  the  three  chemical  reaction 
models  is  shown.  Figures  20  and  21  present  the  OH  number 
density  distributions  along  the  y  =  8  m  line  and  the 
x=-1.5  m  line,  respectively.  As  expected,  the  TCE  model  of 


FIG.  19.  OH  number  density  (molecule/m3)  contours  with  the  rate  of  Xie  et 
al.  (TCE)  at  120  km  altitude  for  the  freestream  velocity  of  5  km/s.  Area 
shown  is  7.4  X  7.5  m. 
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FIG.  20.  Comparison  of  OH  number  density  (molecule/m3)  at  y  =  8  m  for 
three  chemical  reaction  models  at  120  km  altitude  for  the  freestream  veloc¬ 
ity  of  5  km/s. 


case  (1)  predicts  higher  OH  production  than  the  MD/QCT 
model.  With  the  case  (1)  reaction  rate,  the  OH  production 
was  much  more  than  the  other  two  cases  (2)  and  (3)  because 
the  O  +  HCl  reaction  probability  was  significantly  higher,  as 
shown  in  Figs.  2  and  6.  With  the  case  (2)  reaction  rate,  the 
maximum  OH  number  density  is  reduced  by  more  than  50% 
whereas  with  the  MD/QCT  probability,  it  is  slightly  lower 
than  for  case  (2).  There  are  two  factors  that  caused  the 
change  in  OH  production  between  case  (2)  and  (3).  The  first 
factor  is  the  larger  O  +  HCl  collision  cross  section  obtained 
by  the  MD/QCT  method  that  leads  to  more  OH  molecules 
being  produced.  The  second  factor  is  the  lower  MD/QCT 


FIG.  21.  Comparison  of  OH  number  density  (rnolcculc/m ')  at  x=— 1.5  m 
for  three  chemical  reaction  models  at  120  km  altitude  for  the  freestream 
velocity  of  5  km/s. 
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FIG.  22.  Reaction  probability  distribution  for  O+HCl— >OH+Cl  using  the 
TCE  model  (the  rate  of  Xie  et  al.)  for  an  altitude  of  120  km  and  a  freestream 
velocity  5  km/s. 


reaction  probability,  however,  the  effect  of  the  change  of 
collision  cross  sections  was  found  to  be  smaller  than  the 
change  of  reaction  probability.  Thus,  these  two  factors  result 
in  a  slightly  lower  OH  number  density  compared  to  case  (2). 

In  the  shock  region,  the  reaction  probabilities  changed 
dramatically  between  cases  (1)  and  (2).  For  case  (1),  the 
reaction  rate  results  in  TCE  reaction  probabilities  greater 
than  one,  as  cases  found  in  our  previous  work.  (In  the  DSMC 
calculation,  they  are  set  to  one.)  Due  to  the  jet-atmosphere 
interaction,  high  translational  energies  result  in  high  reaction 
probabilities  for  the  TCE  model.  About  10%  of  the  O+HCl 


Reaction  probability 

FIG.  23.  Reaction  probability  distribution  for  O+HCl— > OH + Cl  using  the 
MD/QCT  model  for  an  altitude  of  120  km  and  a  freestream  velocity  of 
5  km/s. 
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collisions  results  in  the  probability  greater  than  one  with  the 
rate  (1)  using  the  TCE  model  as  shown  in  Fig.  3.  In  contrast, 
with  the  rate  of  Xie  et  al.,  case  (2),  most  of  the  O  +  HCl 
collisions  are  predicted  to  have  a  TCE  reaction  probability 
smaller  than  0.4  (see  Fig.  22).  Finally,  the  reaction  probabil¬ 
ity  with  the  MD/QCT,  case  (3),  is  still  lower  than  the  other 
two,  and  more  than  80%  of  the  reaction  had  the  probability 
of  lower  than  0.1  as  shown  in  Fig.  23. 

Using  the  same  VHS  (Bird  values)  collision  cross  sec¬ 
tion  for  the  TCE  and  MD/QCT  reaction  models,  the 
MD/QCT  predicts  higher  reaction  probability  than  case  (2) 
as  shown  in  Fig.  6.  In  contrast,  if  both  the  reaction  and  col¬ 
lision  cross  sections  are  obtained  by  the  MD/QCT  method, 
the  MD/QCT  reaction  probability  becomes  lower  than  the 
TCE  with  the  rate  (2).  At  5  krn/s,  the  MD/QCT  and  the  case 
(2)  TCE  reaction  probabilities  are  very  close,  but  at  3  km/ s 
or  8  km/s,  the  case  (2)  reaction  probability  is  higher  than 
those  obtained  with  the  MD/QCT. 

VI.  CONCLUSION 

Chemically  reacting  flows  of  the  interaction  of  a  jet  with 
a  transitional-to-rarefied  atmosphere  have  been  simulated  us¬ 
ing  the  DSMC  method.  The  objective  of  this  work  was  to 
assess  the  sensitivity  of  the  flow  modeling  to  the  fidelity  of 
chemical  reaction  models  in  hypervelocity  collisions  in 
DSMC.  The  0(3H)  +  HCl(1X+)^0H(2n)  +  Cl(2H)  chemical 
reaction  was  selected  for  study  since  it  is  an  important  reac¬ 
tion  path  in  atmospheric-jet  interaction  flows.  The  MD/QCT 
calculations  were  performed  for  the  O  +  HCl  reaction  using 
the  ab  initio  3A"  potential  energy  surface  of  Ramachandran 
and  Peterson  (RP)  (Ref.  8)  for  both  reaction  and  collision 
cross  sections.  The  approach  of  Tokumasu  was  utilized  to 
calculate  the  viscosity  cross  sections,  and  the  VHS- 
equivalent  collision  cross  sections  were  obtained  with  the 
assumption  of  isotropic  scattering.  Finally,  the  MD/QCT  re¬ 
action  probabilities  were  compared  with  those  obtained  by 
the  TCE  model  with  two  rates. 

These  results  show  that  using  the  same  VHS  collision 
cross  section,  the  MD/QCT  reaction  probability  was  pre¬ 
dicted  to  be  lower  than  the  TCE  probability  with  the  rate  of 
Mahmud  et  al.1  used  in  Ref.  3(1)  and  higher  than  the  TCE 
probability  with  the  rate  of  Xie  et  al.10  (2).  The  TCE  model 
with  rate  (1)  predicted  reaction  probabilities  higher  than  one 
in  the  energy  range  of  interest.  The  VHS  cross  section  was 
found  to  be  very  sensitive  to  the  viscosity  parameters,  obvi¬ 
ously.  Therefore,  to  be  accurate,  the  MD/QCT  collision  cross 
sections  using  the  approach  that  Tokumasu  used  were  calcu¬ 
lated  for  the  range  of  relative  velocities  and  internal  energies 
of  interest. 

The  energy  transfer  between  the  translational  and  inter¬ 
nal  modes  was  also  investigated.  For  the  initial  HC1  internal 
energy  of  5  X  1CT20  J,  translational  energy  was  transferred  to 
the  HC1  internal  modes.  The  average  of  energy  shift,  how¬ 
ever,  was  found  to  be  almost  zero  for  an  initial  HC1  internal 
energy  of  2  X  10-19  J.  The  magnitude  of  the  inelastic  cross 
sections  is  small  compared  to  the  total  cross  section.  The 
dependence  on  the  initial  internal  energy  is  so  small  that  the 
MD/QCT  VHS  equivalent  collision  cross  sections  were  con¬ 


verted  to  fit  the  VHS  with  «=0.39  and  <iref=3.9  at  1000  K. 
These  values  were  used  for  the  MD/QCT  reaction  probabili¬ 
ties,  and  these  probabilities  were  tabulated  and  used  in 
DSMC. 

The  DSMC  simulations  were  performed  for  the 
atmosphere -jet  interaction  flows  for  an  altitude  of  120  km,  a 
freestream  velocity  5  km/s,  and  three  O  +  HCl  chemical  re¬ 
action  models.  The  chemical  reaction  models  did  not  affect 
the  overall  flow  field,  but  affected  the  properties  of  the  prod¬ 
uct  species,  OH.  Although  the  TCE  with  rate  (1)  resulted  in  a 
reaction  probability  higher  than  one,  the  TCE  with  rate  (2) 
and  MD/QCT  models  showed  that  the  O  +  HCl  reaction  has  a 
maximum  probability  less  than  0.4.  Also,  the  OH  production 
was  slightly  lower  with  the  MD/QCT  than  the  TCE  rate  (2) 
model. 

For  low  enthalpy  reactions,  in  hypervelocity  collisions, 
the  TCE  model  with  accurate  Arrhenius  rates  appears  to 
agree  well  with  rigorous  MD/QCT  calculations.  Good  agree¬ 
ment  between  TCE  and  MD/QCT  is  observed  because  this 
reaction  does  not  show  strong  favoring.  However,  in  the  re¬ 
gime  where  the  TCE  reaction  probability  approaches  unity, 
one  should  verify  the  total  collision  cross  section. 

Note  added  in  proof.  A  very  recent  work  by  Camden  and 
Schatz"  provides  important  new  direct  dynamics  results  on 
the  reaction  cross  sections  for  O+HCl. 
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